arXiv:1507.03500v3 [cond-mat.quant-gas] 21 Apr 2016 


Extended Bose-Hubbard Models with Ultracold Magnetic Atoms 


S. Baler,^ M. J. Mark,^’^ D. Fetter,^ K. Aikawa,^’[^L. Chomaz,^’^ Z. Cai,^ M. Baranov,^ R Zoller,^’^ and F. Ferlaino^’^ 

^Institut fur Experimentalphysik, Universitdt Innsbruck, Technikerstrafie 25, 6020 Innsbruck, Austria 
^Institut fUr Quantenoptik und Quanteninformation, Osterreichische Akademie der Wissenschaften, 6020 Innsbruck, Austria 
^Institut fur Theoretische Physik, Universitdt Innsbruck, Technikerstrafie 2IA, 6020 Innsbruck, Austria 

(Dated: April 22, 2016) 

The Hubbard model underlies our understanding of strongly correlated materials. While its standard form 
only comprises interaction between particles at the same lattice site, its extension to encompass long-range 
interaction, which activates terms acting between different sites, is predicted to profoundly alter the quantum 
behavior of the system. We realize the extended Bose-Hubbard model for an ultracold gas of strongly magnetic 
erbium atoms in a three-dimensional optical lattice. Controlling the orientation of the atomic dipoles, we reveal 
the anisotropic character of the onsite interaction and hopping dynamics, and their influence on the superfluid- 
to-Mott insulator quantum phase transition. Moreover, we observe nearest-neighbor interaction, which is a 
genuine consequence of the long-range nature of dipolar interactions. Our results lay the groundwork for future 
studies of novel exotic many-body quantum phases. 


PACS numbers: 67.85.Hj, 37.10.De, 51.60.+a, 05.30.Rt 

Dipolar interactions, reflecting the forces between a pair of 
magnetic or electric dipoles, account for many physically and 
biologically significant phenomena. These range from novel 
phases appearing at low temperatures in quantum many-body 
systems IHIEI, liquid crystals and ferrofiuids in soft condensed 
matter physics Ellll, to the mechanism underlying protein 
folding O. The distinguishing feature of dipole-dipole in¬ 
teractions (DDI) is their long-range and anisotropic charac¬ 
ter a pair of dipoles oriented in parallel will repel each 
other, while the interaction between two head to tail dipoles 
will be attractive. While remarkable progress has been made 
with gases of polar molecules (71 and Rydberg ensembles (H 
comprising electric dipoles, it is the recent experimental ad¬ 
vances in creating quantum degenerate gases of bosonic and 
fermionic magnetic atoms, including Cr IMD and the Lan¬ 
thanides Er CD and Dy ca , which have now opened the door 
to a study of magnetic dipolar interactions, and their unique 
role in Hubbard dynamics of a quantum lattice gas. 

Ultracold Lanthanide atoms with their open electronic f- 
shells, and their anisotropic interactions are characterized by 
unconventional low energy scattering properties, including the 
proliferation of Feshbach resonances (141 . This complexity of 
Lanthanides manifests itself in quantum many-body dynam¬ 
ics: by preparing quantum degenerate Lanthanide gases in op¬ 
tical lattices we realize extended Hubbard models for bosonic 
and fermionic atoms. Here, in addition to the familiar sin¬ 
gle particle tunneling and isotropic onsite interactions (as for 
contact interactions in Alkali) dipolar interactions give rise to 
anisotropic onsite and nearest-neighbor (offsite) interactions 
(NNI), and density-assisted tunneling (DAT) na. Such ex¬ 
tended Hubbard models have been studied extensively in the¬ 
oretical condensed matter physics and quantum material sci¬ 
ence caini, and it is the competition between these uncon¬ 
ventional Hubbard interactions, which underlies the predic¬ 
tion of exotic quantum phases such as supersolids, stripe and 
checkerboard phases |[T8H23]| . 

Here we report a first observation of the unique manifes¬ 


tations of magnetic dipolar interactions in extended Hubbard 
dynamics. These observations are enabled by preparing an ul¬ 
tracold sample of bosonic Er atoms in an three-dimensional 
(3D) optical lattice. It is the control of the optical lattice via 
laser parameters in combination with a fiexible alignment of 
the magnetic dipoles in an external magnetic field, which al¬ 
lows us to reveal and explore the anisotropic onsite and offsite 
interactions. Measurements of the excitation spectrum in the 
Mott insulator state, and of the superfiuid-to-Mott-insulator 
(SF-to-MI) quantum phase transition are employed as a tool 
to detect these interactions and their competitions. 

In our experiment an ultracold dipolar gas of ^^^Er atoms 
is prepared in a 3D optical lattice. The atoms are spin- 
polarized in their lowest Zeeman sublevel ifTTl and feature 
a magnetic moment jl of 1 Bohr magneton. The experi¬ 
ment starts by adiabatically loading a Bose-Einstein conden¬ 
sate (BEC) of about 1.5 x 10^ atoms from an optical dipole 
trap (ODT) into the optical lattice. The lattice is created by 
two retrorefiected 532-nm laser beams, defining the horizon¬ 
tal xy-plane, and one 1064-nm beam, nearly collinear with 
the vertical (z) direction given by gravity (Eig. lA) (Supple¬ 
mentary Materials). The lattice has a cuboid unit cell with 
lattice constants dx^ = 266 nm and d^ = 532 nm, which cor¬ 
respond for Er to the recoil energies Er = Er 3; = hx 4.2 kHz 
and Er^^ =hx 1.05 kHz, h being Planck’s constant. In addi¬ 
tion, the lattice can be controlled by independently changing 
the depths associated with the lattice beams in each direction, 
{sx^Sy^Sz), measured here in units of the corresponding recoil 
energies. The dipole orientation, quantified by the polar an¬ 
gles 6 and (j) (inset Eig. lA), is varied by changing the direc¬ 
tion of the polarizing magnetic field (Supplementary Materi¬ 
als). By changing the lattice depths we can prepare the Er 
atoms in a Mott insulator state, driving a SE-to-MI phase tran¬ 
sition, as described below. 

The dynamics of Er atoms in the optical lattice is described 
by an extended Bose-Hubbard (eBH) model with Hamilto- 
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FIG. 1. (color online) Magnetic dipoles in a 3D optical lattice. (A) Schematic of our lattice geometry, where the lattice constants are indicated. 
The dipole orientation, given by the polarizing magnetic field B, is quantified by the polar angles 0 and 0 with respect to our coordinate system. 

(B) Illustration of the contributing terms in the eBH model: Tunneling matrix element Jij, DAT matrix elements A/^y, onsite interaction t/, and 
NNI Vij. (C) Illustration of the definition of the onsite aspect ratio. (D to F) Calculated values of the DDI-dependent terms as a function of 0 
for (j) = 0° and typical experimental parameters {sx,Sy,Sz) = (15,15,5^) with set by the AR for the cases AR = 1 (red) and AR = 2 (green). 

(C) shows C/^d- (D) gives Vij=x (solid lines) and Vij=y (dotted lines), the NNI for bond direction x and y respectively. (E) shows AJx^dd (solid 
lines) and A/^ (dotted lines), the values of A//^- ^d for hopping direction x and y respectively. The dashed lines indicate the case without 
DDL C/s and Jij are independent on 0 and their values for the two configurations considered are C/s = 3749 Hz (1775 Hz) for AR = 1 (2) and 
Jij = 27 Hz. 
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Here b] (b -) are the bosonic creation (annihilation) operators 
of atoms at site /, nt = bjb- is the associated number operator, 
and (/ j) denotes pairs of adjacent sites. The first term in Eq. 1 
includes the single-particle hopping, with amplitudes Jij re¬ 
flecting the anisotropy of the optical lattice. Interactions mani¬ 
fest themselves in an onsite interactions U, offsite interactions 
Vij (approximated as nearest-neighbor-interaction (NNI)), and 
a density-assisted-tunneling term (DAT) AX , IZTl . All 
terms of the eBH are illustrated in Fig. IB. As discussed in 
the Supplementary Materials (see also ll25]| ). the onsite in¬ 
teraction U and DAT AJtj have contributions from both the 
short-range part of the interatomic interaction (C7s and 
which is proportional to the v-wave scattering length and 
from the long-range DDI (C7dd and dd), which is propor¬ 
tional to JLL^. On the other hand, Vtj originates entirely from 
the long-range DDI. This mechanism for NNI is in marked 
contrast to, for example, Heisenberg spin-spin interaction be¬ 
tween atoms at neighboring sites /, j, which arises from su¬ 
perexchange processes in Hubbard dynamics in second order 


virtual hopping processes ^ jf^/U in the limit of large onsite 
interaction (23 |29l . 

The unique and characteristic feature of our many-body 
system is contained in the angular dependence of U, AJtj, 
and Vij, reflecting the DDI in our eBH model. It reveals it¬ 
self prominently in combination with an anisotropic Wannier 
function at a given lattice site reflecting the local density distri¬ 
bution (Fig. 1C). The aspect ratio AR of the Wannier function 
can be changed by imposing unequal lattice depths {sx^Sy,sf). 
In our experiment Sx = Sy, such that z is the anisotropy axis. 
We deflne AR = where 4 (4 = ly) is the harmonic os¬ 

cillator length along the z (xy) direction of the local atomic 
well oni. The relative weight between the attractive and re¬ 
pulsive contribution to Udd can be tuned by changing the 
dipole orientation relative to the anisotropy axis of the on¬ 
site density distribution, and the AR (Fig. ID). In contrast, the 
NNI Vij is controlled through the orientation of the dipoles 
with respect to the bond direction ij (Fig. IE). Finally, dd 
depends both on the orientation of the dipoles relative to the 
bond and anisotropy axes, and on the AR (Fig. IF). 

We first investigate the impact of the DDI on the onsite in¬ 
teraction (Udd) by performing spectroscopic measurements. 
We prepare our system deep in the MI phase and probe the 
energy gap in the excitation spectrum for different dipole ori¬ 
entations. This energy gap, associated to particle-hole excita¬ 
tions, is U for atoms in singly- or doubly-occupied Mott shells 
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FIG. 2. (color online) Measurement of the onsite interactions. (A and B) Excitation spectrum of the MI state for dipole orientations 0 = 0° 
(A) and 0 = 90° (B). The modulation spectroscopy is performed at {sx^Sy^Sz) = (15,15,52.5), corresponding to AR 1.46 and the remaining 
BEC fraction is measured after ramping down the lattice depths to zero. Erom a double Gaussian fit to the data (solid line) we extract the 
resonant excitation frequency Vex for the U and 2U feature. (C) Vex for the loss feature at U and 2U (inset) as a function of AR for 0 = 0° 
(squares) and 0 = 90° (circles). (D) Difference in the energy gap relative to the two dipole orientations, lAt/^dl. as a function of AR. The error 
bars for all figures are the sum of the SEM and systematic errors (Supplementary Materials). The theoretical model (solid lines) reproduces 
very well the experimental data (C and D) and it also includes the effect of the NNI, which shifts the excitation frequency by up to 3 %. Eor 
completeness, calculations accounting only for the isotropic (contact) interaction are shown (dashed lines). 


and 2U at the border between the two shells (Supplementary 
Materials) ED We excite the MI by applying a sinusoidal 
modulation of frequency Vex on the amplitude of the x-lattice 
beam ll32H34l . When matches U or 2U, we observe a res¬ 
onant depletion of the condensate. We perform the measure¬ 
ment for 0 = 0° and 0 = 90° (Fig. 2A and B) and observe 
that the resonance positions clearly depend on the dipole ori¬ 
entation, consistent with our expectation. 

To further explore this effect, we repeat the measurement 
for different values of the AR (Fig. 2C). For the spherical 
case (AR = 1), we observe that the excitation gap looses 
its angle dependence showing that f/dd averages to zero ll35l . 
As the spatial distribution is deformed towards larger AR, 
we find a clear deviation from the purely contact-interaction 
case (dashed lines), with a smaller energy gap for dipoles at 
0 = 0°, and a larger one for dipoles at 0 = 90°. Our measure¬ 
ment shows that Udd plays a fundamental role in the stability 
of the MI phase: it can either protect the MI phase for the 
dominantly repulsive DDI (0 = 90°) or make it more suscep¬ 
tible to excitations for the dominantly attractive case (0 = 0°). 
The energy difference between the two dipole configurations 
|A(/dd| is shown in Fig. 2D. The observed angle dependence is 
well described within our eBH model (Fig. 2, C and D, solid 
lines), with the only fit parameter. The derived value for 
(2s = 137(1) ao, with ao being the Bohr radius, is consistent 
with previous measurements based on thermalization experi¬ 
ments CEl. 

In principle, the energy gap in the MI phase also depends 


on the NNI between atoms occupying adjacent lattice sites. 
However, in the above described measurement its influence 
is veiled by the much larger Udd- To demonstrate the pres¬ 
ence of the NNI, we design a dedicated measurement scheme, 
which allows to isolate it from the other terms of the eBH. 
The measurement is based on modulation spectroscopy in the 
2D short-spacing lattice plane (xy-plane), where the NNI is 
stronger (Fig. 3A). The key idea is the following: For a sys¬ 
tem with only onsite interactions the energy gap associated 
with the particle-hole excitation does not depend on the direc¬ 
tion of excitation, i.e. on the direction of the modulated beam. 
In contrast, a system including anisotropic NNI will exhibit a 
modification of the energy gap according to the excitation di¬ 
rection as the energy gap equals U — Vij for excitations along 
the bond direction i j. Hence the difference between the two 
resonance frequencies measured by modulating and de¬ 
noted as AVnni/^, directly reveals the existence of the NNI 
as the onsite contribution cancels. Our scheme is illustrated 
in Fig. 3, A to C, for the case 0 = 90°, = 90°. Here, one 

bond of attractive (repulsive) NNI with energy (y^®P) gets 
destroyed during the excitation along (perpendicular to) the 
dipole orientation such that AVnni = Our mea¬ 

surement for two dipole orientations in the plane with (j) = 90° 
(0°) give AVnniA = +74(10)Hz (-87(14)Hz). Remark¬ 
ably, IAVnniI is similar for both values of (j) as expected from 
the symmetry between these two configurations and is close 
to the theoretical expectation hx 91 Hz, as shown in Fig. 3D. 
This set of measurement provides the first observation of the 
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FIG. 3. (color online) Nearest neighbor interactions. (A) Initial system in the MI regime with dipole orientation 6 = 90°, 0 = 90°. Driving 
excitations along y (B) or v (C) leads to two different particle-hole energy gaps U — and U — Here, = Vij with in-plane head- 
to-tail dipole orientations and = Vij for the in-plane side-by-side orientation. The difference between the two resonant energies AVnni, 
here equals to —+ reveals the NNI. (D) Histogram of AVnni for two dipole orientations 0 = 0° and 0 = 90°. The solid lines in 
front of the histograms are the normal distributions of the corresponding data. On the bottom plane the two dashed lines show the theoretical 
expectation values, while the solid lines show the corresponding measured value with the shaded areas indicating the SEM. 


NNI predicted by the eBH model. 

Finally, we study the effect of the anisotropic DDI on 
the many-body phase transition from a SF to a MI phase 
(Fig. 4A). This phase transition is a result of the competition 
between interactions and tunneling (361, therefore also reveal¬ 
ing the influence of the DAT ca. We probe the atomic inter¬ 
ference patterns of the expanding cloud after a sudden release 
from the 3D lattice for different dipole orientations (0=0° 
and 90°). For increasing lattice depths the system enters the 
MI phase, as shown by the disappearance of the interference 
pattern and the increase of the full width at half maximum 
(FWHM) of the central interference peak (Fig. 4, B to D). 
We clearly observe a shift of the phase transition as a function 
of the lattice depth depending on 0 and AR. For large ARs 
(Fig. 4B), the MI phase is favored for 6 = 90° with respect 
to 0 = 0° as expected from the angle dependence of t/dd- For 
the spherical case (AR =1) (Fig.4D), contrary to the naive 
expectation, we also observe the shift of the phase transition, 
which is now inverted, with the MI phase favored for 0 = 0°, 
as compared to the case of large AR. Instead, we find that the 
shift of the phase transition vanishes at AR ^1.2 (Fig.4C). 
This behavior is a direct consequence of the action of the 
anisotropic DAT term dd in the eBH. For a more quantita¬ 
tive analysis, we systematically study the phase transition as a 
function of AR for 0 = 0° and 0 = 90° . We extract the crit¬ 
ical value of the lattice depth, Sc{0), which defines the onset 
of the SF-to-MI phase transition (Supplementary Materials). 
As shown in Fig. 4E, the difference Asc = *s'c(0°) — Vc(90°) 
decreases when lowering the AR, crosses zero at AR ^ 1.2, 
and eventually become negative at even smaller AR. This be¬ 
havior is very well reproduced by our calculations using the 
full eBH model within a mean-field approximation (Supple¬ 
mentary Materials), highlighting the importance of DAT in 
the many-body system dynamics. 


Quantum degenerate gases of magnetic Lanthanide atoms 
in optical lattices offer a new avenue to access the physics of 
strongly correlated systems for both bosonic and fermionic 
Hubbard dynamics in the presence of dipolar interactions, 
while building on the well-developed toolbox to prepare ul¬ 
tracold dense samples, and manipulate and measure these 
atomic gases. We have realized the extended Bose-Hubbard 
Hamiltonian with anisotropic onsite and offsite interactions, 
which reveal themselves in the excitation spectrum and in 
the many-body dynamics of the system. Our results show 
how to control the Hamiltonian terms with the dipole orien¬ 
tation and accomplish the long-awaited observation of NNI 
in Hubbard dynamics. An outstanding challenge for future 
experiments is the preparation of exotic quantum phases for 
bosons and fermions: an example is provided by stripe phases 
due to NNI, which become accessible with present interac¬ 
tion strengths at temperatures in the few nK regime (Supple¬ 
mentary Materials). Dipolar interactions can be increased by 
working with Feshbach molecules of magnetic Lanthanides, 
essentially doubling the magnetic dipole moment ITtI . These 
opportunities offered by Lanthanides to access the multitude 
of many-body phases predicted for dipolar quantum matter are 
complemented by the remarkable experimental developments 
with heteronuclear molecules and Rydberg atoms m 
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the Lise-Meitner program of the FWF. The Innsbruck theory 
group is supported by the SFB FoQuS, by the ERC Synergy 
Grant UQUAM, and by the EU FET Proactive Initiative SIQS. 

















5 


E 

S 


E 

S 


X 

5 

LL 


E 

S 


300 

250 

200 

150 

100 

50 

300 

250 

200 

150 

100 

50 

300 

250 

200 

150 

100 

50 



FIG. 4. (color online) Superfluid-to-Mott-insulator transition. (A) Time-off-flight absorption images of the atomic cloud taken 27 ms after a 
sudden release from the 3D lattice with Sxj^z = ^ during ramp up 5 = 0, = 10, 5 = 22 (A, 1 to 3) and during ramp down ^ = 4, ^ = 40 (A, 
4 to 5). (B to D) The width of the central interference peak is plotted as a function of the lattice depth in the xj-plane for AR = 2,1.28,1, 
with dipoles oriented along 0 = 0° (squares) and 6 = 90° (circles). We extract the phase transition point for each orientation and AR via 
a double-line fit and the visibility (Supplementary Materials). (E) shows A^c = ^c(0°) — *^0(90°) as a function of the AR. The data from the 
visibility (diamonds) and the double-line fit (triangles) shows similar results. The dashed red line is the weighted mean of the two methods. 
The solid black line is the theoretical calculation from a mean held approximation (Supplementary Materials). The dotted line shows the 
expectation without any DDL 


SUPPLEMENTARY MATERIALS 
BEC production 

We create a BEC of about 1.5 x 10^ atoms by 

means of evaporative cooling in a crossed optical dipole trap 
(ODT) dJl. The cloud has typically a BEC fraction above 
80%, which is extracted by a two-dimensional bimodal fit to 
an absorption image of the atomic cloud after a time-of-fiight 
(TOE) of 27 ms Ca. The cloud temperature is estimated to be 
about 70 nK. The ODT is operated at 1064nm and is created 
by two beams, one propagating horizontally and one verti¬ 
cally. The beams cross at their respective focal points. The 
elliptic horizontal beam has a vertical (horizontal) waist of 
about 18/im (117/im) and the elliptic vertical beam has a 
waist of about SSjim (110/im) along (perpendicular to) the 
axis of the horizontal beam. The measured trap frequencies 
are (c0x,C0y,C0^) = 271 x (29.0(6),22.2(4), 165.2(5))Hz. We 
observe a lifetime of the trapped cloud of about 10 s. 

The atomic cloud is spin-polarized in the lowest Zeeman 
sublevel (/ = 6,m/ = —6), where J denotes the total angu¬ 


lar momentum quantum number and mj is its projection along 
the quantization axis. The spin polarization already occurs in 
the magneto-optical trap |[3^ and is maintained in the ODT by 
applying a bias magnetic field with a fixed value of 0.40(1) G. 
As discussed below, the magnitude of this field is kept con¬ 
stant for all the experiments, whereas its orientation is varied 
to set the desired dipole orientation. 

3D lattice setup 

We describe the 3D lattice setup in the coordinate system 
given by the two horizontal lattice beams denoting the x and 
y-axis and the direction of gravity giving the z-axis (inset, Eig. 
lA). The horizontal lattice beams are created by two retrore- 
fiected beams with a waist of about 160 jim and a wavelength 
Ax = Ay = 532 nm. The vertical lattice beam has a waist of 
about 300/im and a wavelength A^ = 1064nm. The result¬ 
ing 3D optical lattice is given by y(x,y,z) = VxCOS^{kxx) + 
VyCOS^{kyy) -\-VzCOS^{kzz), where Vt is the lattice depth in the 
/-direction and ki = 271/Xi the corresponding lattice wavevec- 
tor with i = (x,y,z). Because of the different wavelengths, 
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the atoms experience different recoil energies Er / in the 
plane with respect to the vertical direction. The recoil energies 
given by / {2mXf) are Er^j^ = Er^3; = /z x 4.2kHz 

and Er ^2 = hx 1.05 kHz. Here, h is the Planck constant and m 
the mass of the Er atom. Eor convenience, we give the lattice 
depth in units of the corresponding recoil energy Si = Vi /Er^/. 
The maximum lattice depth we can achieve is {sx^Sy^s^) = 
(30,30,220). Because of the Gaussian profile of the lattice 
beams the atoms experience an additional harmonic confine¬ 
ment. At a typical 3D lattice depth of Sy,s^) = (20,20,20) 
we measure (co^, cOy, co^) = Inx (34(1),31(1),43(l))Hz. 

We note that the vertical lattice beam is tilted from the ver¬ 
tical axis hy 6 = 10(2)° and has an azimuthal angle of (j) = 
5(5)°. This has two consequences: (a) The lattice spacings 
and are modified to d^ = 270(2) nm and d^ = 540(4) nm 
with respect to the A /2 case and (b) the tilt of the wavefront 
of the vertical lattice beam gives rise to an additional poten¬ 
tial difference between neighboring lattice sites along x of 
200(40) Hz due to gravity. While (b) only leads to a broad¬ 
ening of the excitation resonances in the modulation spec¬ 
troscopy measurements, (a) could in principle change the val¬ 
ues of the eBH terms. Therefore, we recalculate them consid¬ 
ering our effective lattice spacings for a typical experimen¬ 
tal condition of {sx^Sy^s^) = (15,15,15). We find that the 
isotropic terms are reduced by 3 % while the anisotropic terms 
can differ between 2-6 %, depending on the dipole orientation 
and the direction of the observed process (see Table 1). This 
gives rise to a downshift of the phase transition point of 
about 1 % for both 0 = 0° and 6 = 90°. However, all these 
shifts are not resolvable within our statistical errors and can 
therefore safely be neglected. 

TABLE I. Difference of the eBH terms between the A/2-spacing 
and the actual spacing given in percentage of the A/2-case for three 
dipole orientations (6 = 90°, (j) = 0°), (6 = 90°, (j) = 0°), and 

e = 0°. 



e = 90°, 0 = 0° 1 0 = 90°, 0 = 0° 1 0 = 0° 

Us 

3% 

Jij=x,Jij=z 

3% 

■^ij=y 

0% 

Um 

6% 

-2% 

2% 

u 

3% 

-2% 

3% 

AVnni 

3% 

2% 

- 

^Jij=x 

3% 

2% 

3% 

^Jij=y 

4% 

2% 

3% 

AJij=z 

3% 

2% 

3% 


Lattice depth calibration and onsite aspect ratio 

To calibrate the depths of the horizontal lattice beams we 
use the standard Kapitza-Dirac diffraction method 1391 . Eor 
the vertical lattice we use the technique of parametric heating, 
in which the atoms are excited from the first to the third lat¬ 
tice band UniED. With these methods, we extract the lattice 


depths with an uncertainty of up to 4%. 

The onsite AR is defined in terms of a Gaussian approxi¬ 
mation to the corresponding Wannier function: AR = Iz/lxj, 

where Ixj = dxj /and lz = dz /are the harmonic 
oscillator lengths associated to the lattice beams along x,y and 
z respectively (Note that we use in our measurements). 

The uncertainty of the AR results from the uncertainty of the 
lattice depths and is about 1 %. 

Because of the non ^-state character of Er atoms in their 
electronic ground state, the atomic polarizability of Er has a 
tensorial contribution, which is about 3 % of the scalar one for 
an off-resonant trapping light l4^ . In our system this effect 
gives rise to a different lattice depth depending on the dipole 
orientation. We carefully studied this effect by calibrating 
each lattice beam for both orientations, dipoles aligned paral¬ 
lel or orthogonal to the lattice beam. Our measurements reveal 
that a parallel orientation gives an up to 4 % deeper confine¬ 
ment compared to the orthogonal orientation. Eor simplicity, 
we account for this effect only by a systematic error in the AR, 
leading for instance to the asymmetric error bars in Eig. 2(C 
and D), and Eig. 4E. 

Loading of the 3D lattice 

Eor our experiments the atoms are adiabatically loaded to 
the 3D lattice by an exponential ramp to the final value within 
150 ms, during which the vertical ODT is linearly lowered 
to zero. To perform modulation spectroscopy, the horizon¬ 
tal ODT is switched off within 1ms after the loading. Eor 
the measurement of the BEG depletion, we exactly reverse 
the described process. In the MI phase we estimate a central 
density of two atoms per lattice site. The external harmonic 
confinement leads to a density distribution with a central dou¬ 
bly occupied Mott shell, consisting of up to 40 % of the to¬ 
tal atoms, surrounded by a singly occupied shell. The exter¬ 
nal harmonic confinement is given by the sum of the ODT 
potentials and the Gaussian profiles of the lattice beams dur¬ 
ing the lattice loading. Eor our typical lattice depth condition 
(sx^Sy^Sz) = (20,20,20), the lifetime of the atomic sample in 
the lattice is 5(1) s. In addition, we observe a heating, which 
leads to a full depletion of the recovered BEG for a holding 
time in the lattice of about 1 s. The origin of this heating is not 
fully understood and might be due to frequency fiuctuations 
of the 532 nm laser source. 

Control of the dipole orientation 

The dipole orientation follows the direction of the magnetic 
field, which we control using three pairs of independent coils 
oriented perpendicular to each other. Each pair of coils is in¬ 
dependently calibrated by performing radio-frequency spec¬ 
troscopy, where resonant excitations to higher Zeeman sub- 
levels can be used as a measure of the actual magnetic field 
at the position of the atoms. The dipole orientation can be 
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Holding time (ms) 

FIG. 5. (color online) Measurement of U by the collapse-and-revival 
technique. The FWHM of the central peak of the interference pattern 
is monitored as a function of the holding time after a sudden quench 
from the SF to the MI phase for an initial dipole orientation of 0 = 0° 
(squares) and 0 = 90° (circles). The latter measurement is vertically 
offset by 300 fim for a better visualization of the two data sets. Each 
point is obtained by two to four independent measurements and the 
shaded region indicates the SEM. The solid lines are the fits of a 
damped sine to the data, used to extract the oscillation frequency and 
hence U. 

changed from 0 = 0° to 0 = 90° and for any value of (j). 
Noise of the ambient magnetic field leads to fluctuations of 
the absolute angles 0 and (j) by about 1° around their set val¬ 
ues. During the evaporative cooling sequence the dipoles are 
aligned at 0 = 0°. Before loading the atoms into the 3D lat¬ 
tice, the dipole orientation is changed to the desired value in 
38ms, while the magnetic-fleld magnitude is kept constant. 
After the release of the atoms from the trap, the magnetic 
fleld is rotated towards the imaging direction (0 = 90° and 
(j) = 160°) to perform standard absorption imaging IT^ . 

Modulation spectroscopy in the MI 

To probe the excitation gap in the MI we use a modulation 
spectroscopy technique ||32j |33l. We sinusoidally modulate 
the power of one horizontal lattice beam with a typical to¬ 
tal amplitude between 30% and 40%, and a modulation time 
between 50ms and 100ms. With this method, we resonantly 
create particle-hole excitations in the system El- These exci¬ 
tations manifest themselves as a resonant depletion of the re¬ 
covered BEC because of the extra energy stored in the system. 
We record the remaining BEC fraction after ramping down 
the lattice as a function of the modulation frequency. The 
resulting loss spectrum is then fitted with a double-Gaussian 
function, whose centers give the excitation frequencies. The 
typical EWHM of the resonant loss features is 1 kHz for ex¬ 
citations using the x-lattice beam and 0.8 kHz for the y-lattice 
beam. The width is mainly determined by the external har¬ 
monic confinement. We note that the difference in width be¬ 
tween the two excitation directions is due to the tilt of the 
vertical lattice beam as discussed above. 

We also measure the onsite interaction by using an al¬ 
ternative method, known as the collapse-and-revival tech¬ 


nique 113. Here, we first prepare the system at the onset of 
the SE-to-MI transition with a lattice depth of (sjc,Sy,s^) = 
(10,10,10) and we then suddenly quench the system to 
(sjc,Sy,s^) = (20,20,40) within 5/is. As a result of the 
quench the system oscillates between the MI and the SE phase. 
Eigure SI shows the evolution of the EWHM of the cen¬ 
tral interference peak as a function of the holding time af¬ 
ter the quench for two different dipole orientations 0 = 0° 
and 0 = 90°. We observe up to four collapses and revivals 
and extract the onsite interaction from the oscillation fre¬ 
quency. Eor 0 =0° (0 = 90°) we measure a frequency 
of 2.07(16) kHz (2.98(5) kHz), which are consistent with the 
value of 2.15(3) kHz (2.77(3) kHz) obtained with the modu¬ 
lation spectroscopy technique. 


Analysis of the NNI 

To derive the NNI we perform a differential measurement 
based on modulation spectroscopy, in which the orientation 
of dipoles is fixed but the direction of excitation is changed 
between the horizontal lattice axes x and y. To explain the 
amount of energy needed to drive a particle-hole excitation 
we consider the situation where the dipoles are aligned with 
angles 0 = 90° and (j) = 90°, as also illustrated in Eig. 3 
(see main manuscript). Here we denote (y^®P) the attrac¬ 
tive (repulsive) value of Vij for the bond direction y (x). At 
the starting configuration (Eig. 3A) the total energy is E’a = 
12 yatt _|_ i2Y^^P^ For an excitation along the y-axis the final en¬ 
ergy of this configuration reads SiS Eb = f/ + 1+ 12y^®P, 
while for an x-excitation itisEc = U +1 ly^^P. From 

this consideration it becomes clear that the difference in en¬ 
ergy E’b—E’c = —y^^^ + y^®P = AEnni purely reveals the NNI. 

Analogously the same consideration can be applied for an 
initial dipole orientation of 0 = 90° and 0=0° leading to 
AEnni = — + y^^^ From the theory we expect fh = 

31.5Hz, y^«//z = -59.5Hz and thus IAEnniI A = 91Hz. In¬ 
cluding the corrections arising from the modification of the 
lattice spacings due to the tilt of the vertical lattice beam (see 
above) |AyNNi|/^ changes to 89 Hz, even closer to our mea¬ 
sured values. 

In Fig. S2A we show two excitation spectra obtained us¬ 
ing the method described above described. The difference be¬ 
tween the centers of the Gaussian fits to the data is found to 
be 72(30) Hz and corresponds to one data point of Fig. S2B, 
where all taken measurements are summarized. We be¬ 
lieve that the fluctuation of AEnni along the data sets is 
mainly caused by relative drifts of the lattice depths dur¬ 
ing a differential measurement. We carefully check for 
systematic errors on AEnni using different initial lattice 
depths or atom numbers, but do not And an effect within 
our measurement resolution. The used lattice depths are 
(sjc.Sy.s,) = (15,15,30), (14,18,30), and (20,20,40). The 
different depths can slightly modify AEnni by maximum 2 % 
which is not resolvable within our error bar. 
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Modulation frequency (kHz) Dataset number 

FIG. 6. (color online) Measurements for the NNL (A) Excitation spectrum with modulation along x (diamonds) and _y (triangles) with 6 = 90° 
and (j) = 90°, forming one differential measurement. Each point is the average of about 5 independent measurements and the shaded region 
indicates the SEM. The solid lines are weighted Gaussian fits to the data. The dashed lines indicate the obtained resonance frequencies. (B) 
Set of differential measurements as presented in (A) for dipoles aligned with 6 = 90° and (j) = 90° (squares), 6 = 90° and (j) = 0° (circles). 
The solid lines are weighted fits with the shaded region being the SEM. The dashed lines are the theoretical expectations. 


Analysis of the SF-to-MI transition point 

The critical value of the SF-to-MI transition, Vc, depends 
on the ratio of the total onsite interaction to the total tunnel¬ 
ing rate. In presence of dipole-dipole interaction (DDI), both 
terms depend on the dipole orientations, imprinting an angle 
dependence on the phase transition point Sc{0), defined by the 
value of the horizontal lattice depth = 5-;^ = ‘^ 3 ;- We study 
the phase transition for 0 = 0° and 6 = 90° and extract the 
difference in the critical point Asc = *s'c(0°) — Vc(90°). In par¬ 
ticular, we ramp up simultaneously the three lattice beams in 
150ms while lowering the vertical ODT to zero. We then sud¬ 
denly switch off all beams, let the cloud expand for 27 ms and 
perform standard absorption imaging. We repeat this cycle for 
various final lattice depths j and various values of the AR. 

To extract we use two methods. The first method (a) 
analyses the increase of the FWHM of the central interference 
peak as a function of s^j (Figure S3, A to C). In general we 
observe a smooth transition from the SF phase (small FWHM) 
to the MI phase (large FWHM) as expected for a trapped sys¬ 
tem with a spatially depending density. A weighted fitting 
function consisting of two smoothly connected lines is used 
to extract the critical depth. The region for fitting starts at 
Vjcj = 5 and goes up to a maximum FWHM of 250 /im for all 
data. The point where the two lines are crossing is interpreted 
as the phase transition point 5'c. We carefully check the de¬ 
pendence of the chosen boundaries of the fit region and do not 
find a significant influence on the qualitative behavior of Avc . 
However, quantitatively Avc can vary by up to ±0.2. 

The second method (b) analyses the visibility. Figure S3, 
D to E, shows the extracted visibility data for the same data 
as in (a). The visibility is calculated from a two-dimensional 
fit consisting of seven gaussians for the central and first- 
order interference peaks and one broad gaussian for the in¬ 


coherent background to the obtained interference pattern (see 
Fig. S2G). The visibility is defined as ^ =A/{A-\-B), where 
A stands for the mean amplitude of the first-order interference 
peaks and B for the mean value of the incoherent background 
computed at the positions of the interference peaks. We ex¬ 
tract y as a function of and fit the whole dataset by the 
phenomenological function = C/{1 ^ exp{a{sxj— 

Sc))) — % (adapted from ll44l i. Here C, a, Sc, and % are fit¬ 
ting parameters, where Vc corresponds to the phase transition 
point. 

For both methods, at large ARs, we observe a clear shift 
of toward higher values for 0=0° compared to 0 = 90° 
(Fig. S3, A and D). This shift gets reduced when lowering 
the AR, vanishes around AR ^1.2 (Fig. S3, B and E), and 
changes sign for AR = 1 (Fig. S3, C and F). 

Extended Bose-Hubbard model from microscopic Hamiltonian 

Here we present the details of derivation of the eBH model 
Eq. I together with the expressions for all its coefficients in 
terms of microscopic parameters of the system (see, e. g. 1451 ). 

The microscopic Hamiltonian of the consider system of po¬ 
larized (magnetic) dipolar atoms has the form: 

Aot=^0±^int, (2) 

where the first term 

Ho = l dr^'\r)[-^+V{r)]'¥{r) (3) 

is the single-particle Hamiltonian with 'F(r) being a boson 
field operator, which describes the motion of an atom with the 
mass m in the optical-lattice potential y(r) = VxCOS^{kxx) ± 
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Lattice depth xy-plane 


FIG. 7. (color online) Derivation of the phase transition point (A to F) The FWHM of the central interference peak and the visibility are 
displayed as a function of the final lattice depths for AR = 2 (A and D), 1.28 (B and E), and 1 (C and F), for dipole orientations 0 = 0° 
(squares) and 6 = 90° (circles). The solid lines show the corresponding fitting functions (see text). Each point is obtained from about 10 
independent measurements. (G) Example of density distribution after TOE (top) and the corresponding 2D fit (bottom) used for the derivation 
of the visibility. 


VyCos^{kyy) -^VzCos^{kzz), and the second term 

dr Jdr'^V"^ (r)'I'''' {r')U{r - r')'P(r')'r(r) (4) 

corresponds to the interatomic interaction. In the considered 
case, the interaction contains a short-range part, which can 
be modeled by a contact potential with the 5'-wave scattering 
length ( 2 s, and the DDI (see, e.g. 1461 ) 



f/(r-r') 


m 


5(r-r') + 


1 — 3C0S^ 0r-r' 

4;r |r —r'l^ 


with 0r-r' being the angle between the relative position of two 
dipoles r — r' and their polarization. 

The Hamiltonian of Eq. 3 determines the single-particle 
band structure. 


2m 


+ V(r)]Mctp(r) = ea(p)Map(r), 


where Map(r) is the Bloch wavefunction corresponding to the 
band a and quasimomentum p from the Brillouin zone (BZ), 
defined hy —n jdi < pijh < nidi with dt ^njki being the lat¬ 
tice spacing along the /-direction, pi the corresponding com¬ 
ponent of p, and £a{v) is the corresponding energy. For our 
purposes it is more convenient to work with Wannier func¬ 
tions = EpeBzexp[-/p(r-R,)]t/ap(r)’ which are lo¬ 

calized at different sites Rt of the lattice and orthogonal to 
each other with respect to both the lattice position i and the 
band index a, / Jr0*^(r)0yj 3 (r) = Using these func¬ 

tions as a single-particle basis in the bosonic field operator, 
where bt^a are the bosonic annihila¬ 
tion operators for particles on the site / in the band a, we 
can rewrite the initial Hamiltonian ^ in terms of the opera¬ 
tors bi^a and bj To obtain the eBH model, we keep terms 


with operators for the lowest energy band only. Note that this 
approximation is legitimate because the interatomic interac¬ 
tion in our case is order of magnitude less than the band gap 
such that the admixture of the higher bands can be neglected. 
From the remaining terms we then neglect those which con¬ 
tain square and higher power of exponentially small spatial 
overlaps of the Wannier functions from different sites (see, 
ll45ll for details). Denoting the operators and the Wannier 
functions for the lowest band as bt, b] and ^/(r), respectively, 
we obtain 

H = - Y^Jij(b]bj + h.c) + Hi (n; - 1) + Vi jHiti j 

{ij) i {ij) 

- ^ AJij[b]bj{ni^nj-l)^h.c] 

(iJ) 


where (/ j) denotes a pair of nearest-neighboring sites. The 
first two terms in this expression correspond to the standard 
Hubbard model with the single-particle hopping amplitude 

Jij = j drip* (r) [- +V ir)](l)j{r) 

and the onsite interaction U = 17^ A- Udd, where comes from 
the contact interaction. 

Us = [drlipii 

m J 


and Uii from the dipole-dipole one, 


Uaa = 


4;r 


- j dr j dr'\^i{r)\- 


1—3 cos^ Gr 


r-r 


./|3 




The third term in Eq. [^corresponds to the NNI with 


V. - 

~ Alt 


- j dr j d r'|<;),(r)|- 


1 — 3cOS^6r_r/ I ^ ,^|2 


r'|3 
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coming from the DDI (the contribution from the contact in¬ 
teraction is proportional to the square of the exponentially 
small overlap and is therefore neglected). Note that the DDI 
also generates interactions Vij beyond the nearest-neighbors, 

which decay as |R; — Ry | . The corresponding terms are ne¬ 

glected in the Hamiltonian of Eq. 5 because they are smaller 
and bring no new qualitative features to the results of the 
present paper. Finally, the fourth term in Eq. 5 describes the 
DAT with the amplitude AJij = s + ^ij,dd resulting from 
the contribution from the contact interaction 

and from the dipole-dipole one 


states with different n (n= I and 2 in our case) are separated 
by the SF phase. 

It should be mentioned that, even though the ME approxi¬ 
mation is known to overestimate the stability of the SF phase, 
here we are interested not in the phase boundary of the SF- 
to-MI transitions itself, but in the relative shift of this bound¬ 
ary when the dipolar polarization is changed from 0 = 0° to 
6 = 90°. In calculating this difference, the MF method is ex¬ 
pected to be much more reliable, as it is demonstrated by a 
good agreement between theoretical and experimental results 
(Fig.4E). 

Observability of the stripe phase 


^ij^dd — ■ 


4;r 


■ j dr'\(l)i{r)\- 


1—3 COS^ 0r-r' 


r-r 


•/|3 




It should be mentioned that in our experiments we also have a 
shallow confining potential Vh(r). It can be taken into account 
by adding the term with 14 ,/ = / <^r|0;(r)pVh(r) to 

the Hamiltonian of Eq. 5. 

The above expressions, together with numerically com¬ 
puted Wannier functions, provide theoretical values for the 
parameters in the eBH model. During the calculations, the 
singularity for |r — r'| ^ 0 in the contributions from the DDI 
was resolved by performing the integration over before 
integrating over |r — r'|. For our experimental conditions, the 
contributions from the DDI is typically few times smaller than 
those from the short-range interaction, but have strong depen¬ 
dence on the form of the Wannier function 0/(r) and on the 
alignment of dipoles relative to the lattice axes (see Fig. 1, D 
to F). 


SF-to-MI transition in the mean-field (MF) approximation 

In the MF approximation (see, e.g. EHU), the ground- 
state wavefunction of the system is written as a product state 
over sites: 


i n=0 


where \n)i denotes the Fock state with n bosons on site i. 

The coefficients are the variational parameters sub- 

2 


jected to the constraint 


C 


{n) 


= 1, which can be de¬ 


termined by minimizing the energy Eg = The 

SF phase is characterized by the local order parameter (bi) = 

{WG\bi\'¥G) = l.n=i ^ 0, which implies that 


are non-zero for several adjacent values of n. In con- 

in) 

trast, in the MI phase C\ are non-zero for only one value of 
n. In a spatially inhomogeneous system (e.g., in the presence 
of a trapping potential 14 ), this value is site-independent, and 
the system has typically a layered structure in which the Mott 


The stripe phase is an example of exotic quantum phases 
induced by the NNI, which is characterized by spontaneous 
translational symmetry breaking along one direction. It can 
be accessed in a deep optical lattice half-filled with atoms, 
when the NNI overwhelms the effects of single-particle tun¬ 
neling and temperature. In this case, the onsite interaction 
is much larger than all the other parameters in the Hamilto¬ 
nian, and prevents two atoms to be on the same lattice site. 
We therefore can consider atoms as hardcore bosons, such 
that the number of atoms on a lattice site can be only zero 
or one. We assume that the dipoles are polarized along the 
v-direction resulting in attractive (repulsive V^^P) NNI for 
the bonds ij in the v(y)-direction, with = 2V. To 

calculate the critical temperature for the stripe phase, we con¬ 
sider typical experimental conditions with = ^' 3 ; = 20 and 
AR =1. In such a lattice, we can ignore the DAT and the 
tunneling in the z-direction, and the single-paticle tunneling 
amplitudes Jtj do not depend on the direction of the hopping, 
Jij = J = hx 20.5Hz. We obtain V = hx 34Hz and the NNI 
value for the bonds in the z-direction is V/S = hx 4.25Hz. 
After neglecting this small coupling in the z-direction, the 
Hamiltonian for each xy-plane can now be written as 

Hxy = -JY,il’]bj + h.c)+Y,{-'2Vnini+i^+Vnini+gy - uni), 

Hi) 

( 6 ) 

where b] and bi are hard-core boson operators and i + Cx 
(i + ey) denotes the neighboring site of site i in the v (y) di¬ 
rection. We also add the chemical potential jl which is chosen 
as /i = —V to satisfy the condition of half-filling (nt) = 1/2. 
To determine the critical temperature of the transition into the 
stripe phases, we perform Quantum Monte Carlo calculations 
based on the worm algorithm for the Hamiltonian of Eq. 6 , 
which is free from the negative sign problem. For the above 
parameters, the calculated value for the critical temperature is 
r^ = 1.4/- 1.5nK. 


* Current adress: Department of Physics, Graduate School of Sci¬ 
ence and Engineering, Tokyo Institute of Technology, Meguro- 
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